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Abstract 

Approximate Bayesian Computation (ABC) enables parameter inference for complex physical systems in 
cases where the true likelihood function is unknown, unavailable, or computationally too expensive. It relies 
on the forward simulation of mock data and comparison between observed and synthetic catalogues. Here we 
present cosmoabc, a Python ABC sampler featuring a Population Monte Carlo variation of the original ABC 
algorithm, which uses an adaptive importance sampling scheme. The code is very flexible and can be easily 
coupled to an external simulator, while allowing to incorporate arbitrary distance and prior functions. As an 
example of practical application, we coupled cosmoabc with the numcosmo library and demonstrate how it can 
be used to estimate posterior probability distributions over cosmological parameters based on measurements 
of galaxy cluster number counts without computing the likelihood function, cosmoabc is published under the 
GPLv3 license on PyPI and GitHub and documentation is available at http: //goo. gl/SmB8EX. 


1. Introduction 

The precision era of cosmology marks the tran¬ 
sition from a data-deprived field to a data-driven sci¬ 
ence on which statistical methods play a central role. 
The ever-increasing data deluge must be tackled with 


new and innovative statistical methods in order to im¬ 
prove our understanding of the key ingredients driv- 


ing our Universe (e.g.. 

Borne, |2009t Ball & Brun- 

ner, 2010; de Souza et al., 20141 

de Souza & Ciardi, 

2015). Given the continuous inflow of new data, one 


does not start an analysis from scratch for every new 
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telescope, but is guided by previous knowledge ac¬ 
cumulated through experience. A new experiment 
provides extra information which needs to be incor¬ 
porated into the larger picture, representing a small 
update on the previous body of knowledge. Such a 
learning process is a canonical scenario to be em¬ 
bedded in a Bayesian framework, which allow us to 
update our degree of belief on a set of model param¬ 
eter^ whenever new and independent data are ac¬ 
quired. 

A standard Bayesian analysis specifies prior dis¬ 
tributions on unknown parameters, defines which pa¬ 
rameter values better describe the relationship be¬ 
tween the model, the prior and the new data, and then 
finds the posterior distribution - either analytically or 
via sampling techniques, as e.g. with Markov Chain 


Monte Carlo (MCMC; [Metropolis et al. [[19531 ). This 
analysis requires a proper construction of the likeli¬ 
hood function, which is not always well known or 
easy to handle. A common solution would be to con¬ 
struct a model for the likelihood (e.g. a Gaussian) 
followed by MCMC, with the expectation that this 
hypothesis is not too far from the true. Nonethe¬ 
less, the challenge of performing parameter infer¬ 
ence from an unknown or intractable likelihood func¬ 
tion is becoming familiar to the modern astronomer. 
Recent efforts to overcome observational selection 


biases in the study of massive ( [Sana et ^|2012| ) and 
not-so massive ( jJanson et al.[|2014l ) stars, to account 
for windowing effects, errors and/or gaps in time- 
series of X-ray emission from active galactic nuclei 
( jUttley et al.[ |2002[ [Shimizu & Mushotzky[ [201 3| > 
and UV emission from stellar coronae ([Kashyap et al. 


2002[ > have been reported. 

The development of ad hoc approaches to this 
problem within astronomy have proceeded indepen¬ 
dently from their long history within the field of pop¬ 
ulation genetics. The latter have ultimately been for¬ 
malized into a rigorous statistical technique known 
as Approximate Bayesian Computation (ABC). The 


*For the purposes of this work, we will only be interested 
in the parameter values of a given model. However, it is im¬ 
portant to stress that in a completely Bayesian approach all the 
elements and hypotheses forming the model can be considered 
part of the prior. In this sense, with the arrival of new essen¬ 
tial information, the Bayesian approach allows for completely 
redefinition of the model itself (Kruschke 201 l[i. 


intuition of ABC dates back to a thought-experiment 
in[Rubin[([1984|), where the basic ABC rejection sam¬ 


pler is used to illustrate Bayes Theorem. Tavare et al. 


( [1997[ ) employs an acceptance-rejection method in 
the context of population genetics, while [Pritchard 


et al.[ ([1999[) presents the first implementation of a 


basic ABC algorithm. Only recently has the ABC ap¬ 
proach been introduced and applied to astronomical 
problems ( [Cameron & Pettitt[ 20 12^ [Schafer & Fr^ 


[man[ [2012[ [Weyant et al.[ [2013^ [Robin et aL| [2014| ). 
This work is part of a larger endeavour natural conse¬ 
quence of such initial efforts. Following the philos¬ 
ophy behind the Cosmostatistics Initiative (COIN)|^ 
we present a tool which enables astronomers to eas¬ 
ily introduce ABC techniques into their daily research. 

The cornerstone of the ABC approach is our ca¬ 
pability of performing quick and reliable computer 
simulations which mimic the observed data in the 
best possible way (this is called forward simulation 
inference). In this context, our task relies on per¬ 
forming a large number of simulations and quanti¬ 
fying the “distance” between the simulated and ob¬ 
served catalogues. The better a parametrization re¬ 
produces the observed data in a simulated context, 
the closer it is to the “true” model. From this simple 
reasoning, many alternatives were developed to op¬ 
timize the parameter space sampling and the defini¬ 
tion of distance. One of such examples is the work of 
[Marjoram et al.[ ( |2003| ), who proposes a merger be¬ 
tween the standard MCMC algorithm and the ABC 
rejection sampling. In astronomy, the method was 
used by [Robin et al.[ ( |20I4| ) to constrain the Milky 
Way thick disk formation. Going one step further, 
[Beaumont et al.[ ( [2009a| ) propose to evolve an initial 
set of parameter values (or particle system) through 
incremental approximations to the true posterior dis¬ 
tribution. The Population Monte Carlo ABC (PMC- 
ABC), method was used to make inferences on rate 
of morphological transformation of galaxies at high 
redshift ( [Cameron & Pettitt| [2012| ) and proved to be 
efficient in tracking the Hubble parameter evolution 
from type la supernova measurements, despite the 


contamination from type II supernova (Weyant et al. 


2013| ). More recently, [Lin & Kilbinger[ ( [2015 ) used 


ABC to predict weak lensing peak counts and [Killedat ' 


I ^http://goo.gl/rQZSAB 
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et al.| ( |2015| > applied a weighted variant of the algo¬ 
rithm to cluster strong lensing cosmology. 

This work introduces cosmoabc]^ the first pub¬ 
licly availablej^Python ABC package for astronom}0 
The package is structured so that the simulation, pri¬ 
ors and distance functions are given as input to the 
main PMC-ABC sampler. In this context, users can 
easily connect the ABC algorithm to their own sim¬ 
ulator and verify the effectiveness of the tool in their 
own astronomical problems. The package also con¬ 
tains exploratory tools which help defining a mean¬ 
ingful distance function and consequently point to 
appropriate choices before the sampler itself is ini¬ 
tiated. 

We first guide the user through a very simple toy 
model, in order to clarify how the algorithm and the 
package work. Subsequently, as an example of cos¬ 
mological application, we show how the machinery 
can be used to define credible intervals over cosmo¬ 
logical parameters based on galaxy clusters catalogues 
Simulations for this example were performed using 
the Numerical Cosmology library (NUMCOSMo;|Vitenti' 


& Penna-Lima[|20T4l ). The connection between cos- 
MOABC and numcosmo is implemented as an indepen¬ 
dent module that can be easily adapted to other cos¬ 
mological probes. 

The outline of this article is as follows. In sec¬ 
tion]^ we give an overview of Bayesian perspective 
and the ABC algorithm, cosmoabc package is pre¬ 
sented in section [^through a simple toy model. Sec¬ 
tion |4] describes in detail how to connect cosmoabc 
and NUMCOSMO to obtain constrains over cosmolog¬ 
ical parameters from galaxy cluster number counts. 
Our final remarks are presented in section]^ 


2. Bayesian approaches to parameter inference 

Statistical inference on unknown parameters is 
often a primary goal of a physical experiment de¬ 
sign. Although it is possible, and at times even de¬ 
sirable, to encounter some unpredictable behaviour 


^https://pypi.python.org/pypi/CosmoABC 
^Shortly after cosmoabc was released 


Akeret et al. 


(2015( 


also presented a Python package for forward modelling through 
PMC-ABC. 

^For similar tools in the context of biology and genetics, see 
e.g. Liepe et al. (2010[i; Oaks (20141. 


among the outcomes of an experiment or a measure¬ 
ment, in many situations the experimentation aims 
at establishing constraints over the parameters of a 
model. In other words, the desire is to use real-world 
data to check prevailing theories. 

In the Bayesian framework, the data are seen as 
the accessible truth regarding a given physical pro¬ 
cess and the model as a representation of our un¬ 
derstanding of such process. This approach is data- 
centred and allows us to update the model whenever 
new information becomes available. In other words, 
our goal is to determine the probability of a model 
given the data. 


piom = 


pm6)pm 

pm 


( 1 ) 


where 6 is the vector of model parameters, D the data 
set, p(9\I)) is called the posterior, the prior, p(9), 
represents our initial expectations towards the model 
and p(D) is a normalization constant. 

In this context, the model parameters themselves 
are considered random variables and each individual 
measurement corresponds to one realization of them. 
Thus, once our prior is confronted with the data, the 
outcome is a posterior probability distribution func¬ 
tion (PDF). Using the posterior distributions we can 
determine credible intervals, which represent our un¬ 
certainty about the model parameter^ For example, 
one may be interested in the most-probable region of 
values for certain parameters. 


2.1. Approximate Bayesian Computation 

The ABC algorithm uses our ability to simulate 
the physical process under investigation to bypass the 
necessity of an unknown or computationally too ex¬ 
pensive likelihood function. It is based on the fol¬ 
lowing crucial elements: 

• a simulator, or forward model, 

• prior probability distributions over the input pa¬ 
rameters p{9). 


®Not to be confused with the frequentist definition confi¬ 
dence interval, where the parameter values are considered fixed 
and therefore, there is no probabilistic interpretation associated 
to them. 
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• a distance function, p(!Di, D 2 ). 

As a simple example, consider the following toy 
model: a given physical process can be probed through 
a catalogue of P observations, £> = {x;, ...,xp}. Our 
model states that this process is driven by a random 
variable, A, following a Gaussian distribution, A ~ 
N(juo,o-o). Thus our goal is to identify credible in¬ 
tervals over po and o-q based on £>. Moreover, our 
prior states that po ^ [P-^P+\ and ctq e [cr_,cr+]. 
Hereafter, we will denote the model parameters as 
e = {p,cr]. 

The main idea behind the ABC algorithm can be 
summarized in three main steps: 

• draw a large number of parameter values, 6\ 
from the prior distribution, p, 

• for each O' generate a simulation, 2)g, and cal¬ 
culate the distance between the observed to the 
simulated catalogues, p' = p{T), Dg), 

• approximate the posterior probability distribu¬ 
tion using the fraction of ff's with smallest as¬ 
sociated distances. 


The above method has been modified and fur¬ 
ther developed in the last decade, generating some 


alternatives to the main algorithm (e.g. Sisson et al. 


2007t IDrovandi & Pettitt||201U [Marin et'ar]|2012t 


Del Moral et al.|20T^ [Ratmann et al.|201^ . One of 


them is presented below. 


2.2. Distance 

In the toy model described above, we can safely 
determine the distance, p, between the measured cat¬ 
alogue D and a simulated one Ds as 


p = abs 




+ abs 


/ cr^, - 0-03 \ 

I /’ 


( 2 ) 


where T> is the mean of all measurements in cata¬ 
logue D and o’© is its standard deviation. Equa- 
tionj^encloses important properties, which should be 
present in any ABC distance function: the distance 
between two identical catalogues is zero and the dis¬ 
tance value increases steeply as parameter values get 
further from the fiducial ones. We emphasis that the 
choice of the distance function is a crucial step in 


the design of the ABC algorithm and the reader must 
check its properties carefully before any ABC imple¬ 
mentation is attempted. 


2.3. Population Monte Carlo ABC 

cosMOABc uses the algorithm proposed by [Beau¬ 
mont et aL] ( [2009a[ ), where successive steps towards 
the posterior are achieved by applying an importance 
(or weighted) sampling in the set of parameter values 
whose distances satisfy a given initial threshold. 

We begin by drawing M values from the prior, 
called particles, {O'} with i e [1, M], such that M » 
N (N is the number of samples needed to character¬ 
ize the prior). For each particle we generate a for¬ 
ward model (simulation) and calculate the distance 
between synthetic and real catalogues p' = pCD, 
From this large set, we keep only the N particles with 
smallest p', which constitute the first particle system 
(St=o) and determine a distance threshold for the next 
iteration (e,=i) as the 75% quantile of all p e St=Q. 
In this initial step, we associate to each particle the 
same weight, = 1.0/A, for 7 g [1, A]. 

In subsequent iterations, t > 0, we perform an 
importance sampling from a popular technique 
where one can draw from a proposal distribution and 
re-weight the particle system so it targets the desired 
posterior distribution. 

The parameter vector resulting from this impor¬ 
tance sampling, 0^^, is used to simulate a catalogue 
and calculate its distance to the observed data, p^y. 
The parameter 0try is stored if ptry < Cf This process 
is repeated until a new set of A parameter values sat¬ 
isfying the distance threshold is completed. For the 
new particle system, the weights are calculated as 


W/ 


p(0i) 




(3) 


where W/ denotes the weight associated to the j - th 
particle in particle system t, p(0j) corresponds to the 
prior probability distribution calculated at Oj , W{_■^ 
is the weight of the i - th particle in particle sys¬ 
tem t - 1 and N(0■’;0■’^_^,Ct-\) represents a Gaussian 
PDlj^ centred in 0 ‘_j, with covariance matrix built 
from St-\ and calculated at Of 


^In general, the Gaussian PDF works well, but can be re¬ 
placed with a different distribution if the parameter space has 
special restrictions, e.g. only takes integer values. 
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Data; D —> observed catalogue. 

Result: ABC-posteriors distributions over the 
model parameters. 

t^O 
K ^ M 

for J = 1,..., M do 

Draw 0, from the prior, p(0). 

Use 0 to generate Ds- 
Calculate distance, p = p(Ds, D). 

Store parameter and distance values, 

Sini ^ {e,p} 

end 

Sort elements in Sini by \p\. 

Keep only the N parameter values with lower 
distance in St=o- 

C,=o <— covariance matrix from St=o 

forL= l,...,Ndo 

I ^ 1/N 

end 

while N/K > A do 
K^O. 
t <— r + 1. 

6t <— 75'*-quantile of distances in . 
while len(»S() < A do 
K ^ K+1 

Draw 00 from St-i with weights Wt-\ . 
Draw 0, from N(0o, C,_i). 

Use 0 to generate 

Calculate distance, p = p(!Ds, D) 

ifp < Ef then 

St^{0,p,K} 

K^O 

end 

end 

for / = 1,..., Ado 

I Wf <— equation (|^. 
end 

Wt <— normalized weights. 

C, <— weighted covariance matrix from 
{S„W,}. 

end 

Algorithm 1: PMC-ABC algorithm implemented 
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Table 1: Glossary for algorithm[T] 


Parameter 

Description 


Observed data set 

2 )s 

Simulated catalogue 

M 

Number of draws for the first iteration 

S 

Particle system 

A 

Number of particles in S 

t 

Time-step (iteration) index 

K 

Number of draws index 

W 

Importance weights 

e 

Distance threshold 

A 

Convergence criterion 

0 

Vector of model parameters 

P(-) 

Prior distribution 

pi.;-) 

Distance function 

Ni*-,0,C) 

Gaussian PDF at • 
with p = 0, cov= C 


Once the new weights are determined, we start 
the construction of a new particle system and the al¬ 
gorithm is repeated until convergence. As pointed 
out by |Beaumont et al.| ( |2009b| ), this is achieved when 
the ABC posterior no longer changes substantially 
with subsequent iterations. Here we consider that the 
system converged when the number of draws neces¬ 
sary to construct a particle system is much larger than 
A (see Algorithmic Table |C and Section |3^ . Each 
iteration brings us closer to the “true” PDF bypassing 
the need of a full likelihood calculation. Moreover, 
as the calculation of one particle is independent from 
the others within each iteration, the algorithm itself 
is more easily parallelizable than a standard MCMC. 

3. COSMOABC 

In COSMOABC, our toy model can be represented 
by a simulation function, 

from scipy.stats import norm 
import numpy as np 

def my_simulation(v): 

””” Toy model simulator ””” 

dist = norm (loc=v [ ^ mean ^ ] , 
scale=v[ ^ std ^ ]) 

11 = dist.rvs(size=v[ ^n^ ]) 


1 

2 

3 

4 

5 

6 

7 

8 

9 


in COSMOABC. 
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11 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 
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return np.atleast_2d(11) . T 

where v is a dictionary of input parameters whose 
keywords mean and std determine the mean and stan¬ 
dard deviation of the underlying Gaussian distribu¬ 
tion, respectively, and n denotes the total number of 
objects in the catalogue. Analogously, a flat prior 
would be written a^ 


from scipy.stats import uniform 

def my_prior(par, func=False): 

’’’’’’Flat prior””” 

gap = par[^pmax^] - par[^pmin^] 
pdf = uniform(loc=par[ ^pmin^ ], 
scale=gap) 
if func == False: 

draw = pdf.rvs() 
return draw 
else : 

return pdf 

with par as a dictionary of input parameters and the 
keys pmin and pmax determining the boundaries of 
the distribution. 

The distance function should be written as 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 
17 


import numpy as np 

def my_distance(d2, p): 

’’’’’’Distance function . ””” 

dmean = np . mean (p [ Mataset 1 ^ ]) - 
np.mean(d2)) 

dstd = np . std (p [ Mataset 1 ^ ]) - 
np.std(d2 ) 


gmean = abs(dmean/ i 

np . mean (p [ Mataset 1 ^ ]) ) 2 

gstd = abs(dstd/ 3 

np . std(p [ Mataset 1 ^ ]) ) 4 

5 

rho = gmean + gstd 6 

7 

return np.atleast_ld(rho) 


and receive as input the simulated catalogue d2 and 
the dictionary p. Notice that the observed catalogue 


is contained in p. So the distance to be calculated is 
between p [' dataset 1' ] and d:^ 

We must store these three functions in one file, 
<f unc_f ile>, and edit the sample input file provided 
within cosMOABC. Each keyword in the sample input 
file is self-explanatory, so here we only emphasis the 
model and prior function parameters 

param_to_fit = mean std 
param_to_sim = mean std n 

mean_prior_par_name = pmin pmax 
mean_prior_par_val = — 2.0 4.0 

std_prior_par_name = pmin pmax 
std_prior_par_val = 0.1 5.0 

mean_lim = -2.0 4.0 
std_lim = 0.1 5.0 

mean = 2.0 
std = 1.0 
n = 1000 

prior_func = my_prior my_prior 

Notice that although the variables mean and std 
are free parameters, we need to provide an initial nu¬ 
merical value, within the constraints allowed by the 
prior. The parameter prior_func stores the prior 
PDF for all the free parameters, in the sequence de¬ 
clared in the variable param_to_f it. Such priors do 
not need to follow the same family of distribution. It 
is possible to define a flat prior for the first parameter 
and a Gaussian one for the second. In that case the 
user input file would include, for example. 


mean_prior_par_name = pmin pmax 
std_prior_par_name = pmean pstd 


prior_func = my_prior gaussian_prior 

considering pmean and pstd as the mean and stan¬ 
dard deviation for the Gaussian prior the second pa¬ 
rameter under investigation. 


^The func argument is needed so we can retrieve a real- —-—^- 

ization and the probability distribution itself. This is used by format was chosen in order to optimize paralleliza- 

cosMOABC in the calculation of the weights. *^*®*'' 
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3.1. Visualizing distance behaviour 

Before we attempt to use the ABC sampler, it is 
important to have an idea of how our distance def¬ 
inition behaves for ditferent combination of model 
parameter values, cosmoabc has a tool which allows 
us to visually inspect the performance of our distance 
definition. The code randomly selects parameter val¬ 
ues from the prior, performs the simulation and cal¬ 
culates the distance for each one of them. These dis¬ 
tances are then plotted as a function of the parameter 
values, one parameter at a time. Ideally, the scatter 
of points in the p x 0 space should present a clear 
minimum in the neighbourhood of the most likely 
parameter value. 

In order to test a personalized distance function, 
do 



Std 

Figure 1: Behaviour of the distance function proposed in our 
toy model as a function of the free parameters mean (top) and 
std (bottom). 


$ test_ABC_distance.py -i <input_file> 

-f <func_file> 


An example of the result of this test for the toy model 
we have been considering is shown in Figure [T| No¬ 
tice that the distance behaves as expected, approach¬ 
ing zero around the fiducial values mean=2.0 and 
std=l. 0 and rapidly increasing as parameter values 
move further away. 

It is worth mentioning that this procedure was 
implemented only to provide the user with an in¬ 
tuition regarding the distance function dependence 
with model parameters. The behaviour illustrated in 
Figure [T] is a necessary but not sufficient character¬ 
istic of an ideal diagnostic. Selecting an appropriate 
distance function is an open and problem dependent 
challenge but it is an active area of statistical research 
(see e.g. [Feamhead & Prangl^ |2012^ [Blum et aL 


2013t [Ratmann et al.[ |2013| ). A deeper investigation 
on the steps leading to an optimal distance definition, 
although very important, is out of the scope of this 
work. 


3.2. Running the ABC sampler 

After we are convinced of the performance of our 
distance function, we can proceed to the ABC sam¬ 
pler run. In cosmoabc, this is done through 


$ run_ABC.py -i <input_file> 
-f <func_file> 


The time necessary for the algorithm to converge 
depends on the efficiency of the simulator, the be¬ 
haviour of the distance function and the number of 
particles in each particle system. We suggest an ini¬ 
tial run with a fairly large convergence threshold, for 
example delta = 0.25. This means that the code 
will run until it is necessary to take 4 times more 
draws than the number of particles in each particle 
system. In order to facilitate debugging and interac¬ 
tion with other codes, for each particle system cos¬ 
moabc outputs parameter values, distance, distance 
threshold, computational time and weights for each 
particle in ASCII tables. 

Once the algorithm converges, it is possible to 
visualize the results with 

$ plot_ABC.py -i <input_file> 

-f <func_file> 

-P T 

This will generate a file containing one snapshot for 
each particle system from t=0 to t=T, as well as plots 
for the evolution of distance threshold, convergence 
criteria and computational time. From this first quick 
test, the user can either be satisfied with the achieved 
result or deeide to continue iterating the sampler. If 
more iterations are required, it is only necessary to 
decrease the parameter delta in the user input file 
and continue from the last eompleted particle system 
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$coiitiiiue_ABC . py -i <input_file> 
—f <func_file> 
-P T 


4. Case study: cosmological parameter inference 
from Sunyaev-Zeldovich surveys 

The current concordance cosmology has been re¬ 
markably successful in explaining the observed prop¬ 
erties of large-scale structures ([Tegmark et al4|2006 


Benson! |2010|>. In this framework, the formation 


of such structures proceeds in a hierarchical man¬ 
ner driven by pressureless cold dark matter, where 
galaxy clusters stand out among the largest bound 
objects observed so far. The development of an un¬ 
derlying theory of cluster formation (see jKravtsov & 
|Borgani[ |2012! for a review), allows us to use the 
abundance of clusters as well as their spatial distri¬ 


bution as powerful cosmological probes (e.g., |All^ 
et al.l|200^ . 


There are, however, a couple of caveats which 
make this an interesting problem for the ABC ap¬ 
proach: the model is not deterministic, in the sense 
that it considers the observed data as a realization of 
a Poisson distribution (analogously to the toy model 
studied before) and the unavoidable modelling of the 
observable uncertainties and errors in both, photo¬ 
metric redshifts and mass estimates (for a mathemat¬ 
ical description we refer the reader to jAppendix A 


Penna-Lima et al.| ( r2014| > and references therein). Us¬ 


ing PMC-ABC surpasses the need to integrate a very 
complex likelihood function and reduces the influ¬ 
ence of initial hypothesis on photometric redshift er¬ 
rors in the estimated posterior PDFs. 

Since there is no previous literature on the appli¬ 
cation of PMC-ABC to this particular problem, it is 
crucial to establish a proof of concept. Thus, here we 
present results from a completely synthetic frame¬ 
work, where the “observed” data, 2), is one instance 
of our forward model. This allows us to provide a 
controlled scenario and to ensure our capability of 
recovering the input parameter values. It also facil¬ 
itates the identification (and quantification) of even¬ 
tual biases in the final ABC-posteriors. 



Figure 2: Distribution of observed features (detection signifi¬ 
cance, and redshift, z) in the “observed” catalogue, D. 


4.]. Simulations or the forward model 

Mock catalogues were generated with the num- 
cosMO librar)!^ which provides a set of tools to per¬ 
form cosmological calculations. The software allows 
a large range of possibilities for input cosmological 
and astrophysical parameters as well as main survey 


specifications (see [Vitenti & Penna-Lima! |2014! for 
a more detailed description). Moreover, it can also 
account for the presence of uncertainties from photo¬ 
metric redshifts and mass-observable relation (here¬ 
after, ^-mass relation, where ^ is the detection signif¬ 
icance) which are crucial for a coherent analysis of 
galaxy cluster number counts. 

Cosmological and astrophysical parameters for 
the fiducial model were chosen in accordance to lRe-l 


ichardt et al.| ( |2013] >: = 0.218, erg 0.807, w = 

-1.01, fib = 0.044, Hq = 71.15 km/s/Mpc, n, = 
0.91, Asz = 6.24, Bsz = 1-33, Csz = 0.83 and Dsz = 
0.24 (see I Appendix A.T] for definitions). Telescope 
characteristics follow the SPT design, with minimum 
and maximum redshifts given by Zmin = 0.3, Zmax = 
1.32, respectively, and survey area AQ = 2500 deg^ 
( jBleem et al.!|2015j ). 

The simulator begins assuming that the total num¬ 
ber of galaxy clusters with z € [z^in, Zmaxl and ^ g 
r^ min , oo) follows a Poisson distribution. It then gen¬ 
erates a realization of this distribution, N^im, and the 
corresponding catalogue {^„Zi}, for i e {1, Ajim} (for 


**'http: //www. nongnu. org/numcosmo/ 
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details in the process see |Penna-Lima et aT||2014[ ap¬ 
pendix B). Here, we investigate the three-dimensional 
space {Qc, os, w} with flat initial priors, Qc € [0.01,0.6], 
cTg e [0.5,1.0] and w e [-3.0,0.0]. All other cosmo¬ 
logical parameters are considered known and fixed at 
the values reported above. 

cosMOABC contains a warp of the numcosmo sim¬ 
ulator which can be accessed through the user input 
file keyword 


1 simulation_func 


numcosmo_sim_cluster 


and an example of the input file with all other options 
tailored for numcosmo simulations is also provided 
within the package. 

Fig. [^displays the static simulated catalogue we 
used as “observed” data in the ^ x z sample space. 
The sample is composed by 671 clusters with z € 
[0.30,1.32] and [5,26]. 

4.2. Distance 

The complexity enclosed in the cosmological sim¬ 
ulations requires some sophistication in designing the 
distance function, cosmoabc has two built-in defini¬ 
tions which proved to be effective in the galaxy clus¬ 
ter counts scenario: quantiles and Gaussian radial ba¬ 
sis function (GRBF) distances. 

The distance_quantile function returns a vec¬ 
tor Q, having L + \ dimensions, where L is the num¬ 
ber of measured feature^ For each feature (column 
in the catalogue), it calculates a few equally spaced 
quantile^ At every quantile, the values of the cu¬ 
mulative distribution functions (CDF) coming from 
simulated and observed catalogues are subtracted and 
the square root of their sum is returned. The last di¬ 
mension accounts for the variability in the total num¬ 
ber of objects. If / is the number of objects in D and 
Is is the number of objects in Ds, the last element of 
p will be 


p_i = max 


abs 1 


, abs 1 


h 

1 


(4) 


of N particles with smaller distances. Once the first 
particle system is constructed, the distance threshold 
€ will also be a L-i-1-dimensional vector. A new set of 
parameter values 0 will only be accepted to populate 
the next particle system if it satisfies the 3 distance 
thresholds independently. 

We emphasis that the this is only a simple and 
computationally fast distance definition which proved 
to be efficient in this synthetic scenario of cosmolog¬ 
ical inference from galaxy cluster number counts for 
the illustrative purposes of this work. Whenever us¬ 
ing ABC in a real data situation, the user must design 
a distance function which preserves these features for 


In the construction of the first particle system, the 
magnitude of this vector, |p|, is used to select the set 


the problem at hand (e.g., see Section 3.3 of |Cameron 
|& Pettitt| ( |20T2| ». 

Figure illustrates the effectiveness of this dis¬ 
tance definition in determining the cosmological pa¬ 
rameters based on SZ flux measurements. The dis¬ 
tance calculations were performed using the cosmoabc 
tool described in section |3.1| however, in order to 
make the visualization lighter, we display binned re¬ 
sults in all three free parameters. In each panel the 
horizontal axis was divided in 500 bins and each dot 
represents the smallest distance found in that bin for 
10"^ draws. From Figure|^we see that the first (com¬ 
parison of CDF over redshift) and second (compar¬ 
ison of CDF over ^) distance elements do present a 
local minimum around the fiducial values for and 
CTg, although the behaviour is much lighter than in 
the previously discussed toy model (Figure [T]). The 
role of the third element (comparison between the 
total number of objects) is to impose an upper limit 
on the free parameter values, since this element in¬ 
creases steadily for ^ 0-28 and erg > 0.86. We 
also see that there is little hope in using this distance 
to constraint w, since there is no significant change 
in behaviour for the three distance elements. 

4.3. Results 

Specific tools are also available for the case of a 
SZ survey using numcosmo. Once all the choices are 
made in the user input file it is possible to run the 
ABC sampler using 


'Tn our case, L = 2, for observed features ^ and redshift. ' ^ run_ABC_NumCosmo . py -i <input_file> 
'^The number of quantiles if defined by the user in the input 
file, through the keyword quantile modes. 
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Figure 3: Behaviour of the quantile distance function in the context of galaxy cluster number counts. Each panel illustrates how 
the elements of the quantile based distance vary as a function of the cosmological parameters for 10"^ random draws from the prior. 
Lines run through parameters and columns through distance elements. 


Analogously, if the user is interested in continuing 
the calculations from a given particle system T, this 
can be done using 


$ continue_ABC_NumCosmo.py 

-i <input_file> 

-P T 

In case a user defined distance or prior is chosen, it 
is necessary to include the -f option followed by the 
name of the function file in both examples above. 
Plots can be generated as shown in section [3^ 

Credible intervals from the ABC-PMC estimated 
posterior distributions are shown in Figure The 
upper panels show 2-dimensions posteriors over the 
free parameters and the bottom panels display the 
profile of each parameter individually. Frames show 
the evolution of the approximated posteriors for con¬ 
secutive particle systems. The first frame merely rep¬ 
resents the initial prior: a flat PDF over all the free 
parameters. The next frame displays results from the 
first particle system (t = 0), where we generate a 
large number of simulations (M = 50000) and kept 


only the 10% with smaller distance. From t = 1 
we clearly see how the posterior evolves and adapts 
through subsequent iterations. The credible intervals 
not only shrink, but also become asymptotically well 
behaved for further particle systems. 

Worth noting, for this example each particle sys¬ 
tem holds N = 5000 particles and the convergence 
was achieved in 9 iterations for delta = 0.01. In 
case a tighter posterior is desirable for and erg, one 
can simply decrease the convergence criteria, letting 
the system evolve for a little longer. If further infor¬ 
mation on w is desired, a more informative distance 
definition should also be usec|^(see Figure [s]). The 
evolution of the distance threshold and convergence 
criteria are shown in Figures]^ andrespectively. 


'^The quantile distance was chosen due to its simplicity and 
low computational cost, cosmoabc also contains a distance def¬ 
inition ( [Appendix B[ ) which accounts for potential correlations 
between two parameters in a catalogue. We advise the user to 
consider the GRBF distance as well as the combination with 
other cosmological probes in case tight intervals over w are de¬ 
sired. 
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Figure 4: Results from coupling cosmoabc to the numcosmo simulator. Frames run from successive iterations of the PMC-ABC 
algorithm. Upper panel: two-dimensional representation of the ABC posteriors in each iteration. Lower left panel: evolution of 
the dark matter density profile. Lower centre panel: evolution of the posterior over erg. Lower right panel: evolution of the PDF 
profile over the dark energy equation of state parameter. 


5. Final remarks 

We presented cosmoabc, a Python implementa¬ 
tion of Population Monte Carlo Approximate Bayesian 
Computation (PMC-ABC) algorithm with adaptive 
importance sampling. Traditional methods of param¬ 
eter inference are useful if the likelihood is avail¬ 
able and feasible to compute. Due to the increasing 
amount of data and their complex modelling in all 
areas of astronomy and cosmology, more and more 
computational power is required in order to explore 
larger parameter spaces whose internal correlations 
can often be impractically complicated or unknown. 
Thus, obtaining a statistical tool which bypasses the 
need of fully evaluating the likelihood is imperative. 
PMC-ABC presents an interesting alternative, cos¬ 


moabc is the first such implementation targeted to the 
astronomy and cosmology community. 

The cosmological simulations are done through 
a connection with the Numerical Cosmology library 
(numcosmo), but the code is flexible enough for user- 
specified distance, simulation and prior functions. 

We stress that ABC is not a substitute for stan¬ 
dard MCMC algorithms when the likelihood is com¬ 
pletely known or easy to calculate. It is a viable alter¬ 
native when we are not able to handle the likelihood 
itself, and thus in situations where a MCMC is not 
feasible. 

In this work, we demonstrated the power of cos¬ 
moabc in estimating posterior probability distributions 
in two situations: a simple toy model and a complex 
cosmological simulation of Sunyaev-Zeldovich sur- 
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Figure 5: Evolution of the distance threshold. The first 
(stars), second (+) and third (x) elements of the quan¬ 
tile distance function were normalized by their respective 
larger values. The horizontal axis runs through all the par¬ 
ticle systems shown in Figure]^ 

vey. In both cases, we demonstrated how cosmoabc 
allows a good approximation of the true posterior 
probability distribution with a fairly simple and user- 
friendly interface. We used a completely synthetic 
environment in order to demonstrate the efficiency 
of the method and to be able to address the accuracy 
of the results. We hope this will be useful not only 
to cosmologists, but to all research areas in astron¬ 
omy where simulations are becoming increasingly 
more accessible and systematics are making likeli¬ 
hood functions even more intractable. 

The code is published under GPLv3 in PyPl and 
Github and documentation can be found in Readthe- 

doc£3 
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Appendix A. Cosmology with galaxy cluster num¬ 
ber counts 


where AQ is the survey area in square degrees and 
H(z) is the Hubble parameter 


The model assumes that the number density of 
halos with mass in the range [M, M+dM] is strongly 
linked to both the matter-energy content of the Uni¬ 
verse and the statistical properties of the initial lin¬ 


ear density contrast field (Press & Schechter, 1974 


Bond et aT||1991t[Sheth & TormenH1999| >. Hence, 
the comoving number density of dark matter halos 
can be written as 


dn(M, z) 
dlnM 


Pm(z) , 1 


do-n 


o-RdlnM' 


(A.l) 


where Pmiz) is the mean matter density at redshift 
z, /(ctr, z) is the multiplicity function ([Tinker et al.[ 


2008), and (t\ is the variance of the linear density 


contrast filtered on the length scale R, 


dk 

—k^P(k,z)\W(k,R)f-, (A.2) 

with a spherical top-hat window function, W(k,R), 
and the linear power spectrum given by 


Hiz) = HolflM + z)^ -t (1 - n„,) (1 + 

(A.7) 

with Qm representing the fraction of total energy den¬ 
sity in the form of matter {Llm), Ho being the Hubble 
constant and w the equation of state of dark energy, 
considered to be constant. 

This framework allows us to relate the mean num¬ 
ber of DM halos within a certain range of mass and 
redshift (equation |A.5| ) to the parameters describing 
the underlying cosmological model (equations jA.lj to 
A.7| ). However, we still need to connect the theoret¬ 
ical redshift z and mass M with their equivalent ob¬ 
servable quantities. We begin by taking into account 
the uncertainties from photometric redshift, Zphoh de¬ 
termination. 

We assume that Zphot follows a Gaussian distri¬ 
bution with mean equal to z and standard deviation 
cr = 0.05(1 -Hz), which we refer to as P(Zphot\z)- Thus, 
the expected number of clusters for a given interval 
of Zphot ^nd M can be written as. 


P(k, z) = Ak^^ T(kf D{zf. (A.3) 


In the expression above the power spectrum depends 
on the spectral index, Ug, the linear growth function, 
D(z), normalized such that D{0) = 1, and on the 


transfer function, T(k) ( jEisenstein & Hu[[1998| >. The 
normalization factor is written as 


A = 


crt 


r 


^k(n,+2)T(kyw\k, 8 ) 


(A.4) 


In order to obtain the mean number of DM halos 
with mass in the range [M, M + dM] and in the red¬ 
shift interval [z, z+dz], we combine the mass function 
with the comoving volume element dV/dz, namely. 


cfN , „ dVdn{M,z) 
-———dzdlnM = -— 
dzd InM dz d\nM 

Assuming a flat universe. 


dzd In M. 


(A.5) 



d^N(M, Zphot, &) 

dZphot^ 


= J" dzPiZphot\z) 


d^N(M,z,e) 
dzd In M 


(A.8) 

where 6 comprises both the cosmological and astro- 
physical parameters (such as those of the mass- ob¬ 
servable relation - see [Appendix A.lj ). 

Estimating the mass enclosed in a given galaxy 
cluster is not a trivial task (see e.g., [Lagana et al.[ 
[2010[ [Giodini et al.[ |2013[ and references therein). 
Traditionally, one requires the recognition of indirect 
signatures carrying such information into observable 
quantities, such as optical and X-ray emissions ( [Birki^ i- 
shaw[[T999t|Carlstrom et al.][2002[). In particular, we 


use a mass-observable relation derived by the SPT 
team, which relies on measurements of the SZ effect. 


Appendix A.l. Cluster mass estimate from Sunyaev- 
Zeldovich effect 

The intracluster medium (ICM) is a hot plasma 
which interacts with photons of the cosmic microwave 
background (CMB) via Compton scattering, causing 
a spectral distortion in the CMB radiation. This is 
known as the SZ effect. The integrated thermal SZ 
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flux, Ysz, is proportional to the total thermal energy 
of the ICM dBarbosa et aH \1996\ |Motl et ahj |2005 ) 
and consequently it is possible to use the SZ distor¬ 
tions on the CMB to estimate the mass of the cluster. 

Due to significant uncertainties in the direct de¬ 
termination of Ysz, we follow here the strategy re¬ 


ported by the SPT (Vanderlinde et al., 2010; Benson 


|et al.[ |2013[ [Reichardt et aH |2013[ ), where an unbi¬ 
ased estimator ^ of the detection significance (sig¬ 
nal to noise ratio) ^ is used as a mass proxy. In this 
context, ^ - 3 and ^ oc YszINint, with 

denoting the noise per resolution element or the inte¬ 
grated noise over several resolution elements for un¬ 
resolved and resolved detections, respectively. More¬ 
over, the adopted mass scaling relation is given by 


f = 4 


SZ 


M 500 


3 X 


Bsz 


E(z) 

£( 0 . 6 ) 


Csz 


(A.9) 


where £(z) = H{z)/Ho, M 500 = (4;r/3) 500pcrit ^ 500 ’ 

with pcrit = SHqISttG as the critical energy density, 
R 500 the radius enclosing 500xpcrit at the cluster red- 
shift, and the scaling relation parameters Asz (^-mass 
normalization), Bsz (^-mass slope) and Csz (^-mass 
redshift evolution) can be determined concomitantly 
with the cosmological parameters. Finally, substitut¬ 
ing the true mass by the unbiased estimator in equa- 
A.8[ the number of clusters with ^ e + d^] 


tion 


and Zphot ^ l.Zphot, Zphot dzphot\ can be expressed as 


cfN(^,Zphot,0) 

dZphotd^ 


-f 


= dzP(Zphot\z) 


f dlnM I dC 


Following [Benson et al.| ( |2013| >, and [Reichardt et al. 
( [2013[ ), we assume 


P(ln ^1 In M)d In f 


1 


X 


exp 


^s/^Dsz 




d^. 


(A.11) 


with Dsz being the log-normal scatter in and 

\2l 


p(md^ 


1 




exp 




d^. 

(A. 12) 


Appendix B. Distance based on Gaussian Radial 
Basis Function 

We describe below a distance defined in terms of 
Gaussian radial basis functions (GRBF), Pgrbf (within 
COSMOABC the function is called distance.GRBF). This 
is not a distance in the mathematical sense, since 
Pgrbf (29, Ds) ^ Pgrbf (^9s , £)), but as the PMC- 
ABC is centred in the observed catalogue, it is enough 
to guide the sampler to the correct posterior distribu¬ 
tion over the evolution of particle systems. 

For a given simulated sample, we compute the 
approximation of its underlying model using a GRBF 
interpolation ( [Tsybakovl [2008[ ), 


Na 


1 / dj c-^ di 

Gr(^, z\D) = 2] - —7== exp 


Inv/MiCj 


(B.l) 


where 


di = 

C = 4cov(2)), 

5 is a scale parameter and cov(£)) is the covariance 
matrix of the observed data. Each element of the sum 


d^N{M,z,e) 
dzd In M 
P(^|^)P(ln^|lnM). 

(A. 10) 


in equation ( [B.l[ ) works as a kernel density distribu¬ 
tion centred in the i-th simulated cluster with (^;, zd- 
Consequently, Gr(^*, z*) is the number of clusters we 
expect to observe having (^*, z*). 

We can re-obtain the total number of objects in 
the catalogue through 

fd?f dzG,((.zlD) ~ ^sim* (B.2) 
Making use of the auxiliary function 


N 


lnf(!Dl!Ds) = ln(GR(zp4|£>s)) - 

7=1 


(B.3) 
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with the index j running through all the data points 
in D dPenna-Lima et al.H2014| >, the distance between 
the “observed” and simulated catalogues is given by 


p(D,Ds) =-2 In 


I fmv) I 


(B.4) 


Therefore, in each iteration, t, we will only accept 
those parameter values whose forward model satisfy 
PtCD, IDs) < €t and the ABC-PMC algorithm can be 
employed normall 3 pj 

The scale parameter 5 (equation |B.2| ) regulates 
our tolerance towards distinct distributions which pro¬ 
duce the same total number of objects in a catalogue. 
Suppose we begin with a large and follow the ABC- 
PMC algorithm reducing €( at each time-step. If s is 
too small the probability of finding parameter values 
satisfying the distance threshold will drop steadily, 
rendering the algorithm unable to further reduce e,. 
On the other hand, if 5' is too large, the density func¬ 
tion will evolve to a very flat behaviour losing all 
information about the underlying distribution of D. 
Thus, 5 must be chosen such that most of the shape 
information in T) is retained, while still being feasi¬ 
ble to reduce 6, until the desired precision is achieved. 
For the specific case outlined here, we found that any 
s e [0.1,0.5] will lead to well constrained and unbi¬ 
ased results. 

In order to call this distance definition from within 
cosMOABC, the input file must include the extra pa¬ 
rameter s and the function definition 


2 distance_func = distance_GRBF 

3 s = 0.15 

4 ... 


The GRBF distance is more time consuming, since 
it takes into account the correlation between the ob¬ 
served features. However, it might be worth to use it 
in highly correlated data scenarios. The current ver¬ 
sion of the NUMCOSMO library includes an ABC sam¬ 
pler using the GRBF distance. 


^^Note that the denominator in equation BA is a scalar inde¬ 
pendent of the mock data. This ensures the distance function 
convergence. 
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